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A SUPERCONVERGENT DISCONTINUOUS GALERKIN 
METHOD FOR VOLTERRA INTEGRO-DIFFERENTIAL 
EQUATIONS, SMOOTH AND NON-SMOOTH KERNELS 

KASSEM MUSTAPHA 



Abstract. We study the numerical solution for Volerra integro-difTcrential 
equations with smooth and non-smooth kernels. We use a h-version discontin- 
uous Galerkin (DG) method and derive nodal error bounds that are explicit 
in the parameters of interest. In the case of non-smooth kernel, it is justified 
that the start-up singularities can be resolved at superconvergence rates by 
using non-uniformly graded meshes. Our theoretical results are numerically 
validated in a sample of test problems. 

Integro-difTerential equation, weakly singular kernel, smooth kernel, DG time- 
stepping, error analysis, variable time steps 

1. Introduction 

In this paper, we study the discontinuous Galerkin (DG) for a nonlocal time 
dependent Volterra integro-difFerential equation of the form 

(1.1) u'{t) + a{t)u{t) + Bu{t) = /(<), 0<t <T with u{0) = uq, 
where B is the Volterra operator: 

(1.2) Bu{t) = I I3{t,s)u{s)ds, 

Jo 

such that, 

(1.3) /3(i, s) = {t- s)"-'^b{s) for aU < s < i < T 

with either a e (0, 1) (weakly singular kernel) or a g No '■— {1, 2, 3, • • • } (smooth 
kernel). Here a, b and / are continuous real valued functions on [0, T]. We assume 
that there exist /u* > such that a{t) > fj,^, for all t £ [0,T]. As a consequence of 
this and the continuity assumptions on the functions a and b ; there exist ^* , /^* > 
such that 

(1.4) ^Ji^ < a{t) < ^* and \b{t)\ < n* for all t e [0, T] . 

For any uq G Ri problem (jl.ip has a unique solution u which is continuously 
differentiable, see for example \X^. However for a € (0, 1), even if the functions a, 
b and / in ()l.ip - ()1.3p are smooth, the second derivative of u is not bounded at 
t = (see [3] and related references therein), and behaves like |u"(t)| < Ct"^^. 
The singular behavior of u near t = may lead to suboptimal convergence rates if 
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we work with quasi-uniform time meshes. To overcome this problem, we employ a 
family of non-uniform meshes, where the time-steps are concentrated near t = 0. 

Various numerical methods had been studied for problem (jl.ip . For instance, 
collocation methods for (jl.ip with a weakly singular kernel were investigated by 
many authors where an O^k'P'^^) {k is the maximum time-step size and p is the 
degree of the approximate solution) global convergence rate had been achieved us- 
ing a non-uniform graded mesh of the form (|2.10p . see for example [U [31 [21] and 
references therein. Spectral methods and the corresponding error analysis were 
provided in [71 [12] assuming that a = 1 and the solution u of (|l.ip is smooth. How- 
ever, for < Of < 1 (that is, the kernel is weakly singular), the spectral collocation 
method were recently studied in [23] where the convergence analysis was carried 
out assuming again that the solution u is smooth. For other numerical tools, refer 
to [231 and references therein. 

In the present paper we shall study the nodal error analysis for the DG time- 
stepping method (with a fixed approximation order) applied to problem (jl.ip . 
Indeed, the DG time-stepping method for (jl.ip when a G (0, 1) has been intro- 
duced in |2], where a uniform optimal ©(fc^"*"^) convergence rate had been shown 
assuming that u is sufficiently regular. In this work, we show that a faster conver- 
gence than 0(fcP+^) is possible at the nodal points. For a weakly singular kernel 
{a £ (0, 1)), we prove that by using non-uniformly refined time-steps, start-up 
singularities near t — can be resolved at 0(/c™'"{P'"+i}+p+1) superconvergence 
rates . Such convergence rates can not be obtained by using the approach given 
in [2]. Very briefly, our proof technique will be carried out in two steps; deriving 
first the global convergence results of the DG method for the dual problem of (jl.ip 
(which is essential for the nodal error but irrelevant for the global error estimates) , 
see Theorem 14.11 Then, we use these results with the orthogonal property of the 
DG scheme for (jl.ip very appropriately (see (jS.ip and Theorem IS.ip to achieve 
nodal superconvergence estimates. For smooth kernels (a € No), we appropriately 
modify our earlier analyses to show nodal superconvergence rates of order 0{k^P^^) 
assuming that the functions a, b and / are sufficiently regular (see Theorem 16. 2p . 

The origins of the DG methods can be traced back to the seventies where they 
had been proposed as variational methods for numerically solving initial- value prob- 
lems and transport problems [TUl [HI [H [SI [5] and the references therein. In the 
eighties, DG time-stepping methods were successfully applied to parabolic prob- 
lems; see for example, [5], where a nodal 0(fc^P+^) superconvergence rate had been 
proved. Subsequently, in [9], a piecewise linear time-stepping DG method had been 
proposed and studied for a parabolic integro-differential equation: 

(1.5) ut + Au + BAu = J in (0,T] x VL with u(0) = v{x) on f7 for a G (0, 1), 

where C M'' is a bounded convex domain, A is a linear self-adjoint, positive- 
definite operator (spatial), with compact inverse, defined in D{A), and where A 
dominates the spatial operator A. A nodal O(fc^) superconvergence rate had been 
derived assuming that b{s) = 1 in (jl.3p . where the error analysis there was based on 
the fact that on each time interval, the DG solution takes its maximum values on 
one of the end points. However, this is not true in the case of DG methods of higher 
order p. The high order time-stepping DG for (jl.Sp was investigated in [15) where a 
global optimal 0(kP~^^) convergence rate had been proved, assuming that the mesh 
is non-uniformly graded. (For other numerical methods for (jl.5p , see [12l [HI [TH [17] 
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and related references therein.) Indeed, our convergence analysis can in principle 
be extended to cover the nodal error estimates from the DG time-stepping method 
of order p, applied to (jl.Sp . 

The outline of the paper is as follows. In Section [2l we introduce the DG 
time-stepping method with a fixed approximation degree p (typically low) on non- 
uniformly refined time-steps with p > I. In Section |3l we give a global formulation 
of the DG scheme, introduce our projection operator, and also provide some tech- 
nical lemmas. In Section 21 we define the dual of the problem (II. ip and then derive 
the error estimates from the discretization by the DG method when a £ (0, 1); see 
Theorem 14. II In Section [SJ we prove our main nodal error bounds. For a € (0, 1), 
an error |C/" — u{tn)\ of order 0(fc"""{p^"+i}+P+i) (i.e., superconvergent of order 

for p = 1 and fcP+2+" foj- p > 2) has been shown provided that the solution u 
of satisfies (j2.7|) and the mesh grading parameter j > {p+ l)/cr; see Theorem 
15.11 In Section m we consider the case a G Nq (in ((LSI)) and thus the kernel is 
smooth. We show a nodal error of order 0{k^''^^) (over a uniform mesh) assuming 
that the solution u of is sufficiently regular, refer to Theorem 16.21 We present 
a series of numerical examples to validate our theoretical results in Section [T] 



2. Discontinuous Galerkin time-stepping 

To describe the DG method, we introduce a (possibly non-uniform) partition of 
the time interval [0,r] given by the points 

(2.1) = tQ <ti < ■■■ <tN =T. 

We set /„ = {tn-i,tn) and fc„ — tn ~ tn^i for 1 < n < iV. The maximum step-size 
is defined as fc = maxi<„<jv kn. We now introduce the discontinuous finite element 
space 

(2.2) Wp = { w : Jjv ^ M : w|/„ ePp, I < n < N} , 

where Jjv = U^^]^/„, and Pp denotes the space of polynomials of degree < p where 
p is a positive integer > 1. We denote the left-hand limit, right-hand limit and 
jump at tn by — v{t~), w" — w(i^) and [w]" = w" — w", respectively. 

The DG approximation U E Wp is now obtained as follows: Given U{t) for 
t S Ij with I < j < n — 1, the approximation [/ G Pp on the next time-step /„ is 
determined by requesting that 



(2.3) 



T rn—l \rn— 1 



U' + a{t)U{t)+BU{t) 



Xdt = UT'Xl-' + / fXdt 
t„-i 



for all test functions X E fp. This time-stepping procedure starts from Ul — uq, 
and after N steps it yields the approximate solution U S Wp for t G J^. 

Remark 2.1. For the piecewise-constant case p = 0, since U'{t) = and U{t) — 
U1 = U'lr^ -. U" for t G /„, the DG method ([O]) amounts to a generalized 
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backward-Euler scheme 
U" - U""i 1 /■*" 

1 pitr 1 piir ^~} ^ pmin{t,tj ) 

= — / fit) dt~— E / - • 

In this case, the nodal and global errors have the same rate of convergence which 
is 0(fc), see [3 Theorem 3.8]. 

For our error analysis, it will be convenient to reformulate the DG scheme (|2.3p 
in terms of the global bilinear form 

(2.4) 

Gn{U, X) = Ulxl+Y^ [[/]" + [U'{t) + a{t)U{t) + BU{t)\ X dt. 

By summing up p.Sp over all the time-steps and using [/£ — uq, the DG method 
can now equivalently be written as: Find U G Wp such that 

(2.5) Gn{U,X) = uoXI+ fXdt VXeWp. 

Jo 

Since the solution u is continuous, it follows that 

Gn{u,X) = uoXI + / fXdt VXeWp. 

Thus, the following Galerkin orthogonality property holds: 

(2.6) Gn{U ~u,X) yxeWp. 

Before stating the regularity property of the solution u of we display in the 

next remark an alternative form of Gn which will be used in our error analysis. 

Remark 2.2. Integration by parts yields the following alternative expression for the 
bilinear form Gn in (|2.4p : 

Gn{U, X) = C/i^ Xi^ - E U1 [X]" 

n=l 



J2 [ [-U{t)X' + a{t)Uit)X + BUit)X] dt. 



Throughout the paper, we assume that the solution u of (|l.ip satisfies: 
(2.7) lu^^^t)] < Gt"-^ for 1 < j < p+ 1 where 1 < (7 < a + 1 

where the constant C depends on j. For instance, if in (jl.ip the function / ~ 
f^^ fi+f^^ f2 for some ni, K2>0 and the functions a, b, fi and /2 are in C^^^[0,T] 
for I < j < p, then (|2.7p holds for cr = 1 + min{Ki, K2,0i}, see [H Section 7.1] for 
more details. 

We notice from (|2.7p that \u^^^t)\ is not bounded near t = for j > 2. Hence, 
to compensate the singular behavior of u near t = 0, we employ a family of non- 
uniform meshes, where the time-steps are concentrated near zero. Thus, we assume 
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that, for a fixed 7 > 1, 

(2.8) fc„ < C^ktl-^/^ and t„ < C^t„_i for 2 < n < TV, 
with 

(2.9) c-yfc'^ < fci < C^k''. 
For instance, one may choose 

(2.10) tn = {n/Nyr for < n < TV. 

Under the assumptions (|2.7p - (|2.9p . we show in Theorem 15.11 that the error |J7" — 
u{tn)\ is of order fc7o-+min{p,i+Q!}^ for 1 < n < A''. So, we have a superconvergence of 
order /jP+i+™in{p.i+"} provided 7 > {p+ l)/a. However, for a quasi-uniform mesh 
(i.e., 7=1) our bound yields a poorer convergence rate of order ;c'^+n"n{p,i+Q}^ 

3. Projection operator and technical lemmas 

In this section we introduce a projection operator that has been used various 
times in the analysis of DG time-stepping methods; see [24], and state some pre- 
liminary results that are needed in our convergence analysis in the forthcoming 
sections. 

For a given function u G C[0,T], we define the interpolant II^u e Hp by 

(3.1) n"it(t^) = u(i„) and / (u - H^u) w = V'yePp_i(/„) 

and for 1 < 71 < iV. From [Tni Lemma 3.2] it follows that 11" is well-defined. 
To state the approximation properties of H", we introduce the notation 

= sup \(j){t)\ for any (j) e C(i„_i,t„). 

Theorem 3.1. There exists a constant C, which depends on p such that: 

(i) For any < q < p and u\j^ G H'^^^{In), there holds 

\n-u - up dt < Ckl^+^ I lu^'+i^p dt for 1 < n < N. 
t„-i Jt„-i 

(ii) For any < q < p and G H'''^^(In) H C(/„), there holds 

||n-u-it||?^ < Cfc29+i /" " |u(«+i)|2dt forl<7i<iV. 

Proof. For the proof of the first bound, we refer to [191 Section 3] or [24l Chapter 
12, Page 214]. For the second bound, see [20l Theorem 3.9 and Corollary 3.10] or 
[21 Equation (12.10)] . □ 

The following two technical lemmas are needed in our derivation of the error 
estimates. The first lemma has been proved in [9j Lemma 6.3]. 



Lemma 3.2. If g e ^2(6, T) and a G (0, 1) then 



pT / pi \ ^ T^a pT pt 

/ / {t-sr-^g{s)ds] dt<— {T-tr-^ / 



g {s)dsdt. 



The next lemma is the following Gronwall inequality; see [9l Lemma 6.4]. 
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Lemma 3.3. Let {aj}jLi and {bj}jLi be sequences of non-negative numbers with 
0<bi<b2<---<bN- Assume that there exists a constant K > such that 

n „t . 

an<bn + K^aj / ' (i„ - t)""^ dt for 1 < n < TV and a G (0, 1). 

Assume further that S = < 1. Then for n = l,---,N,we have an < C6„ 

where C is a constant that solely depends on K , T , a and 5. 

Throughout the rest of the paper, we shah always implicitly assume that the 
maximum step-size k is sufficiently small so that the condition (5 < 1 in Lemma [ 
is satisfied. More precisely, following Lemma lL2l we shall require that 

4T" ( k" <l. 



4. Error analysis of the dual problem 

This section is devoted to deriving error estimates for the DG method applied 
to the dual problem of the Volterra integro-differential equation (11.11) . The main 
results of this section (more precisely. Theorem 14. ip play a crucial role in the proof 
of the superconvergence error estimate in section [5] 

Let z be the solution of the dual problem 

(4.1) - z' -\-a{t)z{t) + B*z{t) =0 forO<t<T, with z(T) = zt, 

where B*v{t) = /3{s, t)v{s) ds {B* is the dual of the integral operator B) . 
Since z has no jumps and since 

T 

[-v{t)z'{t) + a{t)v{t)z{t) + Bv{t) z{t)\ dt 

T 

v{t){-z'{t) + a{t)z{t) + B*z{t)) dt = 0, 

the alternative expression of Gat given in Remark 12.21 vields the identity 

(4.2) Gjv(w, z) ^ v^ZT forah v e G[0, T] . 

(G(0,T] denotes the space of continuous functions on [0,T]). Let Z G Wp denote 
the approximate solution of (14. ip given by 

(4.3) GNiV, Z) = V!^ZT V F e Wp . 
Hence, the following Galerkin orthogonality property holds: 

(4.4) Gn{V,Z ~ z)^Q VFeWp. 

At this stage, the main aim is to estimate the error Z — z in the L2-norm. First 
it is good to notice that (|4.4p is a discrete backward analogue of (|2.6I) . Since it is 
more convenient to deal with a discrete forward problem, we introduce the functions 
z{t) = z{tM — t) and Z{t) — Z{tM — t) and then, (14.41) can be rewritten as; 

(4.5) Gn{Z^~z,V)^{) VFeWp; 

where Gat is defined as in (|2.4p but with a{t) :— a{tM ~ t) in place of a[t) and 
/3(ijv — s,tN — t) in place of f3{t, s). The finite dimensional space Wp is defined as 
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Wp but on the reverse mesh: = to < ti < ■ ■ ■ < In, where ti — ti-i + ki with 

ki = kjy+l-i- 

Setting C = n~z — z and 9 = Z — Il~z where n~ is the interpolant operator 
defined as in (13.11) . but on the reverse mesh. Then (14. 5p imphes that 

(4.6) GNi0,V)^~GN{C,V) yveWp. 

By the construction of the interpolant we have C(i-) = for all n > 1 and hence, 
using the alternative expression for Gn given in Remark [2.2l and J^" ^ V'{t) dt — 
(by definition of the operator n~), 

(4.7) Gjv(c, y) = Y. [mmv{t) + B^mit)] dt 

n=l-^*"-i 

where 

^C(i) = / P{tN-s,tN-t)as)ds. 
Jo 

In the next theorem we estimate the error between z and Z. 

Theorem 4.1. If z is the solution of the backward VIE (j4.ip . and if Z E Wp is 
the approximate solution defined by (|4.3|) . then 

\z-Z\^dt<Ck^"+^\zT\^ 

Jo 

provided that 

(4.8) f " \0{t)\^ dt < G f " \C{t)\^ dt . 
Jo Jo 

Proof. From the decomposition: Z — z = + 9, the triangle inequality, and (|4.8p . 
we have 

r^N piN piN 

(4.9) / \z-Z\^dt= \S^Z\^dt<C \C\^dt. 
Jo Jo Jo 

Thus, the task reduces to bound the right-hand side of (|4.9|) . Starting from the 
relation z{t) = z{tM — t) and recalling that z satisfies (|4.ip . it is clear that z solves 
the VIE: 

i' + a{tN -t)i{t)+ [ (3{tN ~ s,tN -t)S{s)ds = for < i < T, 
Jo 

with z(0) = zt- Hence, an application of (|2.7p for a — a + 1 with z in place of u 
gives 

(4.10) \z'{t)\+t'-''\i"{t)\+t^-''\~z"\t)\<C\zT\. 
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Now, Theorem 13.11 on the reverse mesh (with C in place of 11 u — u) and ()4.10p 

yield 

(4.11) 

E/ m?dt<cY,kt \~z"(t)\^dt<cY,~kii t'"-'\zT\'dt 

n=2"^*"-i n=2 ,i=2 

N N 



n=2 
N 



+2u |2 



n=2 



and on (0, ii), we notice for 1/2 < a < 1 that 

\C{t)\^ dt < C~kt r \i"{t)\''dt<Ck% f\^"~^\zT\^dt<C\zTfk'}+^", 



and for < a < 1/2 that 

FinaUy, combine (|4.9p and (|4.11l) - (|4.12p . we obtain the desired resuh. □ 

In the next lemma we prove the applicability of the assumption (j4.8p . 
Lemma 4.2. For 1 < n < N , we have 

\0{t)\^dt <C \Cit)fdt 

Proof. Choosing V ^ 9 on (0,t„) and zero elsewhere in (|4.6p and (|4.7I) . then using 
the alternative definition of Gat in Remark 12.21 and 9'6 = {d/dt)\9\'^ /2, we observe 
that 

l^(i-)P + l^(?J)P + El[^]'l' + 2 / Ht)\9{trdt 

= -2 / " \d{t)(:{t) + BC(i) + ^{t) dt. 

Jo L J 

So 

/•tn (-tn /-tn 

a(t)|^(t)|'di < / d{t)\9{t)\\C{t)\dt+ / |6C(0 + ^S^WI I^WM*- 



We use the geometric-arithmetic mean inequality \xy\ < ^ — h ^ (valid for any 
e > 0) we find that 



mmwcit)] dt < ^ / ^/W)\o{t)\m\dt 

<\ r miom'dt+fi* I m{'dt 
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and thus 



3 



(4.13) -/ ~a{t)\e{t)\'dt < ^i* \c{t)\'dt+ \Bat) + Be{t)\\e{t)\dt. 

4 JO Ja Jo 

We employ the Cauchy-Schwarz inequahty, again the geometric- arithmetic mean 
inequahty, and Lemma 13.21 fwith T ~ 

in pin I't 

\BC{t)e{t)\dt<^i* / {t- sY-^\C{s)\\e{t)\dsdt 



JO Jo 

<l£- f"a{ty/^\e{t)\ f\t- s)"-^d{s)^/^\C{s)\dsdt 

JQ JO 
' it 2 \ 1/2 / J \ 1/2 



-fcU m\dsj dtj \^J^ ait)\0{trdtj 

< (Q"/" (^j\t-sr-'~a{.sy/'\C{s)\d.s^ dt+\fj ~a{t)m)?dt 



J 7 a{s)\as)\'dsdt+- I dit)\eit)\' dt 



Similarly, we notice that 







\Be{t)0{t)\dt 

<*-^(—) r {In-tY-^ i Ks)\e{s)\'dsdt + -l h{t)\e{t)\' dt. 



Inserting the above bounds in (|4.13p implies that 



2 



d{t)\0{t)\' dt 

, \ 2 n „tj 



<C /*" |C(0pdi + 4^ f^) V r ~a{t)\9{t)\^ dt. 

Jo a Vm*/ ~{Ji,-i Jo 

Therefore, the desired result now immediately follows after applying of the Gronwall 
inequality in Lemma 13.31 and using the assumption (|1.4p on the function 5, (instead 
of a) . □ 



5. SUPERCONVERGENCE RESULTS 

In this section, we study the nodal error analysis of the DG solution U defined 
by (j2.3p with C/*^ ~ uq. We derive error estimate of the DG solution, giving rise 
to superconvergence algebraic rates. Our analysis partially relies on the techniques 
introduced in Chapter 12] for parabolic problems. 

Theorem 5.1. Let a e (0, 1) in il.S]) . Let the solution u of problem (jl.ip satisfy 
the regularity property \2.1\ and let U e Wp be the DG approximate solution defined 
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by p.Sp withp > 1. In addition to the mesh assumption (|2.8p and (|2.9|) . we assume 
that kn > kn-i for I < n < N. Then 

• forp^ 1, 

max |C/"-u(i„)|<Cfcx 1^7' ^ " ^ " ^/'^ 

l<n<A'' " ^ ^' - [P^ 7> 2/ct 

• and for p > 2, we have 

max |[/"-M(t„)| < Cmax{l,logn}fc"+i X P^"' ^ ^ ^ + l)/'^ 

1<"<A'' ~ y^'l- I' B ; \/cP+\ 7>(p+l)/cr. 

Proo/. From (|473l) . (|4?2|) . (fZ!6)) and (|44|) (recall that 77 = Il^ii - w), we observe that 

(C/i^ - u{tN))zT = Gn{U, Z) - Gn{u, z) 

= Gn{u, Z — z) = Gn{ii^ z — Z). 

The alternative expression for Gn given in Remark 12.21 and the equality 77 (i") = 
show that 

(5.2) GNiViZ - Z) = Sin + 52N, 
where 

Sin = - 

To bound Sin and (52Ar, we start from the regularity property (j4.10p and the relation 

z{t) = z{tN — t), and get 

(5.3) \z'{t)\ + {tN - tY"^\z"{t)\ + {tN - t)2-"|z"'(<)| < C\zt\ . 
For p = 1, the orthogonality property of 11^ yields 



^/ Tj{z-Zydt and S2N = iait)r]{t) + Br]{t)) {z - Z){t) dt. 
,=i"'t"-i -'0 



SiN^Y. f" v{t)z\t)dt^Yl /" r7W[^'(i)-2'(i«)]cfi 

= H / / z"{s)dsdt 



and hence, with the help of (|5.3p we have 

|<5iiv| <C||ry||j„^fc„ / |z"(t)|di 
(5.4) »=i 

<Cfc||77l|,7« / {tN-tr-^\zT\dt = CkU\j,\zM/a. 
Jo 

For p> 2, again the orthogonality property of 11^ gives 

^ pt„ N-l t„ 

SiN = -Y. V{t) z'{t) dt=Y, r^{t) [Ii+z'{t) - z'{t)] dt 



ptN pin 

+ / rj{t) / z"{s)dsdt 
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where II+z' is the discontinuous, piecewise-hnear interpolant of z' defined by 

n+z'(0 := z'(t„_i) + \ t - K-i) for t e In 

with := Jj z'{t) dt denote the mean value of z' over the subinterval /„. 

Elementary calculations show that, for t ^ In, the interpolation error has the 
representation 

U+Z'{t)-Z'{t)= f {s^t)z"'{s)ds+^—^ f {tn-sfz"\s)ds, 



and so, by (|5.3p . 

|<^i^| < CM J, E / Wmdt + CkN / \z"{t)\dt 



<cMj.\zT\[J2''n [tN-tr-^dt+kM {tN-tr-ut 



< C||77|U„|zT|(A;"+Mog(iw/fcjv) + 
where in the last step we used; fc„ > /c„_i for n > 1, and 

AT-l 



1 Jt^-l 

< C E / (^w - ty^ dt < Ck^+'^ / (tN - t)'^ dt. 



To bound 52N, we use ()1.4p . integrating, applying the Holder's inequality and then, 
using Theorem 14.11 we notice that 



\52N\<^i*MJ. \\Z{t)-z{t)\ + j^{t-sr-'ds\z{t)-z{t)\j dt 

<M*hiij«(i+t^/«) r \z{t)~z{t)\dt 

Jo 

< CMj. \Zit) - Zit)\'dt^ < Cfc"+lh|l,;JzT| . 

Using Theorem l3.1l the regularity assumption p.7p . and the mesh assumption (|2.9p . 
we get 



\j, <Cki \u'{t)\^ dt <Cki r 2 dt = C-^— < Ck'^", 
Jo Jo 2a — 1 



and for n > 2, we use p.Sp instead of ()2.9p and obtain 



hiii<ce+^ / \u^'+'\t)\'dt 



< Cky+^ r i2-2p-2 < ^^2p+2^2.-2p-2 < ^ ^2p+2^2.-(2p+2)/7 ^ 



n n 

t„-i 



Finally, combine the above estimations from 6in and S2N, and recalling (|5.2p and 
(|5.ip yield the desired bound for n — N. For the nodal error at any time step tng 
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with 1 < 710 < A^, we follow the above steps with no in place of iV, which will then 
complete the proof. □ 

6. Super-convergence analysis for smooth kernels 

In this section, we handle the nodal super-convergence error analysis of the DG 
scheme (12. 3p for problem when a G No (so the kernel is smooth). We use a 
uniform mesh with step-size k = T/N where k is assumed to be sufficiently small. In 
our analysis, we follow the derivations of Sections 0] and [S] with some modifications. 
We assume that the functions a, h and / are sufficiently regular such that the 
solution u of (jl.ip satisfies \u^^\t)\ < C (and consequently |5(-^^(t)| < CIztI) for 
1 ^ J < P + 1 with t £ (0, T]. Thus, from Theorem 13. II we notice that for ti > 1, 



(6.1) Ml<Ck^P+' \u^P+'\t)\^dt<Ck^P+^. 

We start our analysis by deriving the error involved in approximating the solution 



z 



of the backward VIE (liT 



Theorem 6.1. If z is the solution of the backward VIE (|4.ip . and if Z G Wp is 

the approximate solution defined by (|4.3p . then 



/*" \z - Z\^ dt < Ck^P+^\zT\^ . 
Jo 



Proof. First, we recall (|4.13p (over a uniform mesh) 

4 



(6.2) 1/ dit)m\ut<f,* I \at)\^dt+ 1 \BC{t) + B9{t)\m\dt. 

Jo Jo 



Using the fact that a—1 > 0, and the Cauchy-Schwarz and the geometric-arithmetic 
mean inequalities, we observe 

" \BC{t)0{t)\dt < n* l^t°"^\e{t)\ [ \C{s)\dsdt 



"'0 Jo 

<^ f" t''~^&{ty/^\9{t)\ f d{sy/^\c{s)\dsdt 

M* Jo Jo 

t \ 1/2 



< — [ " t°"H{ty/'^\e{t)\ ( [ d{s)\cis)\^ ds) dt 

M* Jo \Jo J 



<{ — ] / t'"-' d{s)\Cis)\'dsdt+- &{t)\9{t)\Ut 
Jo Jo 4 Jo 

f2a / \ 2 ptn -1 ptn 







Similarly, we notice that 



\Be{t)e{t)\dt 



<[^^ a(s)|6l(s)|2dsdt-l-i^ "a(t)|6l(t)|2df. 
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Inserting the above bounds in (|6.2|) yields 

a{t)\d{t)\'^ dt 

<C J^'^ m\'dt + t"^(^^^ pj^'' t"-Ut J\is)\9{s)\' ds. 
Since one can show by induction on a (a G Nq) that 

r t""' dt = -[f^ - = — [r - (j - 1)"] < k'^r-' < k, 

an apphcation of the standard discrete Gronwall Lemma gives 

\9{t)\^ dt < C [ ^ \C,{t)\^ dt iovl<n<N. 
Jo 

Hence, (|4.9p is vahd now and therefore, we obtain the desired resuh after noting 
that 

\m?dt<cY,k^p+^ \z^^p+^\t)\ut<cep+^\zT\'' . 

□ 

In the next theorem we study the nodal error analysis of the DG solution U 
defined by with [7° = itQ. 



Theorem 6.2. Let a G No m hl."^) . Let the solution u of problem he suf- 

ficiently regular and let U G Wp he the DG approximate solution defined by ()2.3p 
with p>l. Then we have 



ma^ |[/!!-u(t„)| < Ck^P+\ 

1<71<N 



Proof. We follow the steps given in the proof of Theorem 15.11 however we use the 
new bounds of 6in and S2N derived below. The orthogonality property of n~ gives 

SiN^-Y^ / V{t) z'{t) dt = Y T]{t) [Uz'{t) - z'{t)] dt 

where Hz' G Wp-i is defined by: for 1 < n < A^, 

Uz' (t-) ^ z' (tn) and [ {z' - Uz') v dt ^ V w G Pp-2(/n) • 
Hence, from Theorem l3.11 there exists a constant C, which depends onp such that: 

lln^' - z'\\j^ < Ck^P-^ f " Iz'^P+^^^dt < Ck^P\zT\^ for 1 < n < 
and thus, using (|6.ip we obtain 

N 

\Sin\ < CkJ2\\v\\i,M^'{t) - z'(t)|U„ < Ck'P+'\zT\ . 
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For the bound of S2N, we use (|1.4|) . integrating, applying the Cauchy-Schwarz 
inequahty and then, using (|6.1|) and Theorem 16.11 we notice that 

\S2N\<^^*\\v\\J. r (\zit)-z{t)\+ f\t-sr-'ds\zit)-z{t)\) dt 



\ ^0 



<^i*MjJl + t%/a) / \Zit)-z{t)\dt 
Jo 

<Ch|U„n \Z{t)-z{t)\^dt\ <Ck'P+'\zT\. 



□ 



7. Numerical examples 



In this section, we present a set of numerical experiments to demonstrate the 
obtained theoretical error estimates and also, to justify the validation of the DG 
scheme (|2.3p for a wider class of intcgro-differential equations. 

Throughout, we consider problem with T — 1, the initial data uq = and 
b{t) = l/r(a) (Here, T denotes the usual gamma function.). Recall that, u denotes 
the exact solution of and U is the DG solution defined by (j2.3|) using 2* (i > 1) 
subintervals, that is = 2* . 

7.1. Example 1. Choosing the coefficient a{t) and the source term f{t) such that 
the solution u of p.ip is given by 

(7.1) u{t)=e+^e-K 

For a G (0, 1), we notice that near t — 0, u" is not bounded, however u is smooth 
away from t = Q. So, we employ a time mesh of the form (j2.10p for various choices 
of the mesh grading parameter 7 > 1 to verify the results of Theorem l5.ll 

Since the exact solution (|7.ip behaves like i""*"! as i 0+, we see that the 
regularity condition ()2.7p holds for a = a + 1. Thus, from Theorem 15.11 and by 
ignoring the logarithmic factor, we expect 

\\U-u\\node-^ max \U'^-u{tn)\ 
l<n<N 

_ fo(fc7(a + l)+min{p,a+l}) for 1 < 7 < (p+ l)/(a+ 1), 

~ |0(fcP+i+™{p,"+i}) for 7 > (p + l)/{a + 1). 
Case 1 Choosing a{t) — 1, thus 

(7.2) f{t) ^{a + l)ec-' + y (-1)4-^^^ + " + ''^ 



i\ r(2 + 2a + i) ' 

i—Q ^ ' 

To illustrate the theoretical results of Theorem 16.21 we choose a — 2 and so, the 
memory term and the solution u are smooth. As expected, the numerical results 
in Table [1] demonstrate nodal errors of order ©(fc^^+i) for p = 1,2, 3. 

In Tables [IHH we displayed the nodal error \\U — u\\node over the mesh (|2.10p 
with iV = 2* and for different values of 7 when a — 0.2 and a — 0.5 (So |u^-'''(i)| is 
not bounded near t = for j > 2.). Results shown in these tables confirm that the 
best convergence rate we can achieve is 0(fcP+i+"""{p^"+i}) and thus our theoret- 
ical results in Theorem 15.11 are sharp in terms of the convergence order. However, 
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i 


P = l 


p-2 


p^3 


2 
3 
4 
5 


3.953e-05 2.622 
5.430e-06 2.864 
7.063e-07 2.943 
8.991e-08 2.974 


1.675e-07 4.934 
5.391e-09 4.958 
1.712e-10 4.976 
5.409e-12 4.985 


1.928e-10 6.949 
1.537e-12 6.971 
1.199e-14 7.002 
1.003e-16 6.902 



Table 1 . The nodal error and the convergence rate over a uniform 



mesh with N — 2' subintervals when a = 2 in ()7.1|) - (|7.2|) . 



i 


7 = 


1 


7 = 1.25 


7 = 1.4 


6 


6.839e-08 


1.886 


3.991e-08 


3.076 


4.883e-08 3.085 


7 


1.522e-08 


2.168 


4.746e-09 


3.071 


5.767e-09 3.082 


8 


3.118e-09 


2.286 


5.668e-10 


3.066 


6.836e-10 3.076 


9 


6.147e-10 


2.342 


6.796e-ll 


3.060 


8.136e-ll 3.070 



Table 2. The nodal error and the rate of convergence for Case 1 



when a — 0.2 and p = 1- 





i 


7 = 1 


7 = 4/3 


7 = (p + 2.5)/3 




6 
7 

8 


5.19e-10 3.08 
6.28e-ll 3.04 
7.73e-12 3.02 


8.23e-12 4.09 
5.01e-13 4.04 
3.10e-14 4.01 


4.43e-12 4.45 
2.00e-13 4.46 
9.01e-15 4.47 


p = 3 


4 
5 
6 


2.36e-09 3.13 
2.83e-10 3.06 
3.47e-ll 3.03 


1.40e-10 4.07 
8.60e-12 4.03 
5.33e013 4.01 


1.80e-ll 5.47 
4.01e-13 5.48 
8.90e-15 5.49 



Table 3. The nodal error and the convergence rate for Case 1 
when a = 0.5 and p = 2 , 3. 



it indicates that in practice we can relax the restriction on the mesh grading expo- 
nent 7. We conjecture that 7 > (p + 1 + min{p, a-|-l})/(cr-|-a-|-l) sufHces to ensure 
0{kP~^^~^'^^'^^P'°'~^^^) convergence. More precisely, we observe 0(fc('^"'"""'"^^''')-rates if 
1^7< (p+l + min{p, a + l})/{a + a + 1). In Tabled! we have chosen a — 0.2 
in (|7.ip - (|7.2p and the DC solution U G Wi (i.e., the approximate solution is a 
piecewise linear polynomial). An 0{k^'^~^"^^'''^) (i.e., 0{k^'^'^)) convergence rate 
has been observed if 1 < 7 < 3/(cr + a + 1) and 0{k^) if 7 > 3/((7 + a + 1). In 
Table H we considered a = 0.5 and U &Wp where p = 2 or 3. An 0(A:('"+"+i)t) 
convergence rate has been demonstrated ifl < 7 < (p-|-2-|-a)/(CT + Q; + l). Finally, 
we chose 7 > (p + 2 + q;)/(o' + a + 1) in Table H and we realized that the order 
of convergence almost matched the one given in the last column of Table [3] where 
7 = (p + 2 + a)/(CT + a + 1) (i.e., the order of convergence did not exceed p + 2 + a 
for p > 2 as the theoretical results suggested). 

Case 2 Choosing a{t) = + 1 and so 

(7.3) f{t) = (« + l)re-* + i2"+ie-* + i^^+i y (-1)^- ^(2 + « + ») 



2—0 
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i 


Error 


Rate 




6 


4.6578C-12 


4.4511 


p= 2, 7 = 5/3 


7 


2.1039e-13 


4.4685 




8 


9.3953e-15 


4.4850 




3 


1.1677e-09 


5.4859 


p = 3, 7 = 2 


4 


2.5328e-ll 


5.5269 




5 


5.6022e-13 


5.4986 



Table 4. Nodal errors and convergence rates for Case 1 when 
a ~ 0.5 and p = 2 , 3. 



i 


7 = 1 


7 = 1.25 


7 = 1.4 


6 
7 
8 
9 


1.633e-07 2.305 
3.205e-08 2.350 
6.220e-09 2.365 
1.205e-09 2.367 


9.644e-08 3.024 
1.184e-08 3.026 
1.454e-09 3.026 
1.787e-10 3.024 


1.233e-07 3.024 
1.514e-08 3.026 
1.858e-09 3.026 
2.284e-10 3.024 



Table 5. The nodal error and the rate of convergence for Case 2 



when a — 0.2 and p — I. 





i 


7=1 


7 = 4/3 


7 = b + 2.5)/3 


p = 2 


6 
7 

8 


1.55e-09 3.01 
1.92e-10 3.01 
2.39e-ll 3.01 


2.62e-ll 4.01 
1.63e-12 4.01 
1.02e-13 4.00 


6.94e-12 4.40 
3.21e-13 4.43 
1.47e-14 4.45 


p = 3 


4 
5 
6 


5.90C-09 2.62 
8.22e-10 2.84 
1.08e-10 2.93 


3.76C-10 3.68 
2.60e-ll 3.85 
1.71e-12 3.93 


1.63e-ll 5.49 
3.59e-13 5.50 
8.07e-15 5.47 



Table 6. The nodal error and the convergence rate for Case 2 



when a — 0.5 and p = 2 , 3.. 



In Tables [5] and |6] we displayed the nodal error \\U — u\\node over the mesh 
(|2.10p with N = 2^ and for different values of 7. Again, we observe convergence of 
order 0(fc('^+"+i)'') ifl<7<(p+l + niin{p,Q! + l})/(a + a + I) and of order 
(3(fcP+i+mm{p,a+i}j -f ^ > {p+l+mm{p, a + 1}) /(o-+a+ 1) for different polynomial 
degrees p. 

7.2. Example 2. In this example we demonstrate that the nodal superconvergence 
results of Theorem 15.11 arc still valid even if a{t) = in (|1.1[) (so the assumption 
(jl.4|) is not satisfied) with a e (0, 1). 

In this case, (jl.ip reduces to the following (scalar evolution or fractional wave 
equation, see [UlIIS]) time-dependent problem: for a e (0, 1), 

/■* (t - sl"-i 

(7.4) u' + / /- — u(s) ds = f(t) for < t < T with u(0) = . 

Jo r(a) 

The piecewise linear {p = 1) DC method for (|7.4p had been studied extensively in 
[T5] . However, for p > 2, the stability and convergence analyses of the DC method 
for (|7.4p are more difficult and it will be a topic of future research. 
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i 


7=1 


7 = 1.25 


7 = 1.5 


6 
7 

8 
9 


9.11e-08 2.33 
1.76e-08 2.37 
3.37e-09 2.39 
6.40e-10 2.39 


1.70e-08 2.79 
2.38e-09 2.84 
3.24e-10 2.87 
4.34e-ll 2.90 


2.59e-08 2.76 
3.67e-09 2.82 
5.05e-10 2.86 
6.83e-ll 2.89 



Table 7. The nodal error and the convergence rates for Example 
2, when a = 0.2 and p = 1. 





i 


7 = 1 


7 = 4/3 


7 = (p + 2.5)/3 


p = 2 


5 
6 
7 


3.94e-09 3.03 
4.89e-10 3.01 
6.09e-ll 3.00 


1.27e-10 4.32 
7.89e-12 4.01 
4.92e-13 4.00 


1.77e-10 4.47 
7.89e-12 4.48 
3.51e-13 4.49 


p = 3 


4 
5 
6 


2.17e-09 3.00 
2.72C-10 3.00 
3.40C-11 3.00 


1.36e-10 4.00 
8.49e-12 4.00 
5.31e013 4.00 


2.23e-ll 5.35 
5.74e-13 5.28 
2.50e-14 5.12 



Table 8. The nodal error and the convergence rate for Example 



2 when a = 0.5 and p = 2 ,3. 



Using the Mittag-Leffler function E^{x) — X]^o^^/r(l m^-Y write 

the exact solution as 

u{t)= f ^„+i(-s"+i)/(t-s)ds. 

JQ 

Choosing a source term /(t) = (a + l)t", we find that 

(7.5) u{t) = -Via + 2)Y: - r(« + 2) (1 - £;.+i(-t"+^)) . 

p—i ^ ^ 

Since the exact solution of (17. 4p behaves like t"''^^ as t — >■ 0+, we see that the 
regularity conditions (|2.7p hold for any a = a + 1. For p = 1 (that is, piecewise 
linear DG method) , the numerical results shown in Table [7] demonstrate a nodal 
superconvergence rate of order 0(fc'''('^+"+^)) for 1 < 7 < 3/((J + a + 1), and of 
order 0{k^) for 7 > 3/(cr + a+l). However, for p > 2, the numerical results 
shown in Table [5] illustrated a nodal error estimates of order 0(fc''(P+^+") (that is, 
0(fc7(<T+P+i)) for 1 < 7 < (p + 2 + a)/((T + a + 1), and almost of order 0(fcP+2+") 
for 7 > (p + 2 + a)/(CT + a + 1). 
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